clear
version 14.0
set more off
set scheme s2mono
capture log close
log using "..\logs\precip results.txt", text replace

* Specify the year of cash rent data to use in the regression
local reg_year 2013

* run the do-files to define programs used in this code
run "..\codeSetup\soils transform - program.do"
run "precip results - program.do"

local bootstrap_reps 1000

*------------------------------------------------------------------------
* Potential soils
*------------------------------------------------------------------------
local soils rootznaws slope om sar gypsum cec7 ec  ksat clay_perc silt_perc ///
			bulkDensity soc0_150 rootznemc ph_less6 ph_greater7dot5 hydgrp* text_*

local soils_nonlin rootznaws slope_* om_* sar_* gypsum_* cec7_* ec_*  ksat_* clay_perc_* silt_perc_* ///
			bulkDensity_* soc0_150_* rootznemc_* ph_less6 ph_greater7dot5 hydgrp* text_*

*------------------------------------------------------------------------
* Controls to include in the regression
*------------------------------------------------------------------------
local deficit ppt
*local surplus water_surplus
local edd dday34C
local soils_controls rootznaws soc0_150_4 bulkDensity_4 ec_2 ph_less6 ph_greater7dot5 slope_1     
*local soils_controls $SelSoils

local knots_`deficit' 4
*local knots_`surplus' 4
local knots_gdd 3
local knots_`edd' 3

*------------------------------------------------------------------------
* Regression
*------------------------------------------------------------------------

use "..\dataAnalysis\rent_climate_data", clear

gen ln_cash_rent=ln(cash_rent_non`reg_year')


ricardian_program_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("AVG") wild(0)

ereturn display
matrix beta_main=e(b)

est restore precip_results

			

*------------------------------------------------------------------
* Graph predicted rent with average soils across precip
*------------------------------------------------------------------
preserve
	foreach var of varlist gdd `edd' `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist gdd `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	qui predict ln_rent_hat, xb
	qui predict ln_rent_hat_se, stdp
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	sort `deficit'
	graph twoway (rarea rent_hat_low rent_hat_hi `deficit', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' `deficit', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `deficit', lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid)  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("Precipitation (inches)") ytitle("Cash Rent ($/acre)") 
	
	graph export "..\figures\precip_predicted_rent_precip.pdf", replace
	
restore


*------------------------------------------------------------------
* Predicted rent across GDD
*------------------------------------------------------------------
preserve
	foreach var of varlist `deficit' `edd' `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist `deficit' `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	qui predict ln_rent_hat, xb
	qui predict ln_rent_hat_se, stdp
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	sort gdd
	graph twoway (rarea rent_hat_low rent_hat_hi gdd, fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' gdd, msymbol(Oh) mcolor(gs11)  msize(small)) ///
	(line rent_hat  gdd, lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("GDD") ytitle("Cash Rent ($/acre)") 
	
	graph export "..\figures\precip_predicted_rent_gdd.pdf", replace
	
restore


*------------------------------------------------------------------
* Predicted rent across EDD
*------------------------------------------------------------------
preserve
	foreach var of varlist `deficit' gdd  `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist `deficit' gdd {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	qui predict ln_rent_hat, xb
	qui predict ln_rent_hat_se, stdp
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	sort `edd'
	graph twoway (rarea rent_hat_low rent_hat_hi `edd', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' `edd', msymbol(Oh) mcolor(gs11)  msize(small)) ///
	(line rent_hat  `edd', lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("EDD") ytitle("Cash Rent ($/acre)") 
	
	graph export "..\figures\precip_predicted_rent_edd.pdf", replace
	
restore

			
simulate _b, seed(2030) reps(`bootstrap_reps') : ricardian_program_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("AVG") wild(1)
save "..\temp\precip_uncertain_reg", replace 
			
mean _all
qui corr _all, cov
matrix var=r(C)

matrix colnames var = rel_total_chg rel_precip_chg rel_heat_chg share_precip share_heat 
matrix rownames var = rel_total_chg rel_precip_chg rel_heat_chg share_precip share_heat


ereturn post beta_main var, dof(14)
ereturn display		






*-----------------------------------------------------------------------------------------------------
* Estimate Damages across all climate models - RCP 4.5
*-----------------------------------------------------------------------------------------------------
local model_list "ACCESS bcc BNUESM CanESM2 CCSM4 CESM1BGC CNRMCM5 CSIROMk3 inmcm4 IPSLLR IPSLMR MIROC5 MIROCESM MIROCCHEM MPIESMLR MPIESMMR MRICGCM3 NorESM1M"

* dataset to store results across climate models (only climate uncertainty)
tempname memhold
tempfile results
postfile `memhold' str15 model rel_total_chg rel_precip_chg rel_heat_chg ///
					share_precip share_heat using `results', replace
* dataset to store results across boostrap reps of all climate models (regression and climate uncertainty)
tempfile bootstrap_results

foreach model of local model_list {
use "..\dataAnalysis\rent_climate_data", clear
gen ln_cash_rent=ln(cash_rent_non`reg_year')
qui ricardian_program_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("`model'") wild(0)

ereturn display
matrix beta=e(b)
post `memhold' ("`model'") (beta[1,1]) (beta[1,2])	(beta[1,3])	(beta[1,4])	(beta[1,5])

simulate _b, seed(2030) reps(`bootstrap_reps') : ricardian_program_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("`model'") wild(1)

capture append using `bootstrap_results'
save `bootstrap_results', replace
			
mean _all
qui corr _all, cov
matrix var=r(C)

matrix colnames var = rel_total_chg rel_precip_chg rel_heat_chg share_precip share_heat 
matrix rownames var = rel_total_chg rel_precip_chg rel_heat_chg share_precip share_heat


ereturn post beta var, dof(14)
ereturn display		

}

postclose `memhold'
use `results', clear
summ, detail
save "..\temp\precip_uncertain_climate", replace

use `bootstrap_results', clear
save "..\temp\precip_uncertain_climate_reg", replace

log close




*-------------------------------------------------------------------
* Create graphs of uncertainty
*-------------------------------------------------------------------
set more off
use "..\temp\precip_uncertain_reg", clear
foreach j in total precip heat {
	summ _b_rel_`j'_chg, detail
	foreach q in 50 75 25 95 5 {
		scalar `j'_reg_p`q'=r(p`q')
	}
}

use "..\temp\precip_uncertain_climate", clear
foreach j in total precip heat {
	summ rel_`j'_chg, detail
	foreach q in 50 75 25 95 5 {
		scalar `j'_climate_p`q'=r(p`q')
	}
}


use "..\temp\precip_uncertain_climate_reg", clear
foreach j in total precip heat {
	summ _b_rel_`j'_chg, detail
	foreach q in 50 75 25 95 5 {
		scalar `j'_both_p`q'=r(p`q')
	}
}


clear
set obs 3
gen uncertainty=_n
foreach q in 50 75 25 95 5  {
	gen total_p`q'=. 
	gen precip_p`q'=. 
	gen heat_p`q'=. 
}


foreach j in total precip heat {

* Regression uncertainty
foreach q in 50 75 25 95 5  {
replace `j'_p`q'=`j'_reg_p`q' if uncertainty==1
}


* Climate uncertainty
foreach q in 50 75 25 95 5  {
replace `j'_p`q'=`j'_climate_p`q' if uncertainty==2
}

* Climate & Regression uncertainty
foreach q in 50 75 25 95 5  {
replace `j'_p`q'=`j'_both_p`q' if uncertainty==3
}

* Format the mean effects to display on the graphs
local `j'_mean_effect=`j'_reg_p50
local `j'_mean_effect : di %9.2f ``j'_mean_effect'
}


* Relative change in total rent
twoway rspike total_p75 total_p95 uncertainty, pstyle(p1) || ///
rspike total_p25 total_p5 uncertainty, pstyle(p1) || ///
rbar total_p50 total_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar total_p50 total_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(-0.3 0.95 "`total_mean_effect'", place(center))

graph export "..\figures\precip_rel_chg_total.pdf", replace

* Relative change due to precip
twoway rspike precip_p75 precip_p95 uncertainty, pstyle(p1) || ///
rspike precip_p25 precip_p5 uncertainty, pstyle(p1) || ///
rbar precip_p50 precip_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar precip_p50 precip_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(-0.05 0.95 "`deficit_mean_effect'", place(center))

graph export "..\figures\precip_rel_chg_precip.pdf", replace

* Relative change due to heat
twoway rspike heat_p75 heat_p95 uncertainty, pstyle(p1) || ///
rspike heat_p25 heat_p5 uncertainty, pstyle(p1) || ///
rbar heat_p50 heat_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar heat_p50 heat_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(-0.2 0.95 "`heat_mean_effect'", place(center))

graph export "..\figures\precip_rel_chg_heat.pdf", replace

